The goal of this document is to examine “escaped muons” that do not decay in the storage region and do not enter a calorimeter. Some of these muons may leave the storage ring altogether.
I generated art files of simulated events using mdc2 with different magnetic fields turned on or off (scenarios). We want to examine,
# Load necessary libraries
library(reticulate) # Access to python
library(stringr) # string functions
library(parallel) # Parallel processing (built-in to R)
library(dplyr) # data analysis
library(purrr) # Functional programming
library(readr) # Read formats
library(ggplot2) # ploting
library(rgl) # 3D plots
library(gridExtra) # Arranging plots
library(testthat) # testing
#
library(readGallery) # Read art Gallery files
# Put all readGallery::useDataProduct calls here
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2ringsim::GeantTrackRecord>')
Samples are located on the Fermilab dCache in the directory /pnfs/GM2/scratch/users/lyon/arr_20170321.
# Function to properly alter paths to accomodate XRootD
# e.g. convert /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
xrootd_ify <- function(aPath) paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(aPath, '/pnfs', ''))
Determine the locations of the data files
# Determine locations of our files
system("ssh lyon@gm2gpvm04.fnal.gov 'ls /pnfs/GM2/scratch/users/lyon/arr_20170321/*/*.root'", intern=T) -> arrFiles
arrFiles %>% xrootd_ify() -> arrXrdFiles
arrXrdFiles
[1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"
[2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"
[3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"
[4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
[5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"
[6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"
[7] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"
[8] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"
[9] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"
[10] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"
[12] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"
Pull out the scenaios (these will be in the same order as the file names)
scenarios <- str_match(arrFiles, 'arr_20170321/.+unified_(.+)_cyl')[,2]
scenarios %>% unique()
[1] "everything" "noDipole" "noInflector" "noKickerQuads" "noKicker" "noQuads"
Let’s look at the contents of these files (they’re all the same, so we’ll just do one) gm2 v7_04_00 will have product_sizes_dumper.
stringToRun <- str_interp(
"ssh lyon@gm2gpvm04.fnal.gov 'source /cvmfs/gm2.opensciencegrid.org/prod7/g-2/setup
setup gm2 v7_04_00 -q prof ;
product_sizes_dumper ${aFile} | grep ::'",
list(aFile = arrFiles[1])
)
#
system(stringToRun)
1410762186 0.456 gm2truth::MagnetIronAndCryoArtRecords_artg4_magnetIronAndCryostat_mdc2arr.
660570363 0.214 gm2truth::TrackingActionArtRecords_artg4__mdc2arr.
637307490 0.206 gm2truth::TrajectoryArtRecords_artg4__mdc2arr.
240738334 0.078 gm2ringsim::GeantTrackRecords_artg4_KeepStepsAction_mdc2arr.
118035827 0.038 gm2truth::RingArtRecords_artg4_Ring_mdc2arr.
12078805 0.004 gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector_mdc2arr.
4482185 0.001 art::RNGsnapshots_randomsaver__mdc2arr.
4086898 0.001 gm2truth::RingTrackingPlaneArtRecords_artg4_RingTrackingPlanes_mdc2arr.
1373178 0.000 gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector_mdc2arr.
1185772 0.000 gm2truth::PhaseSpaceArtRecord_artg4__mdc2arr.
36971 0.000 art::TriggerResults_TriggerResults__mdc2arr.
The tracking action data has birth and death data for every Geant track.
We’ll capture this information for the primary muon and its first generation daughters. See gm2dataproducts/mc/actions/track/TrackingActionArtRecord.hh.
We need a readGallery reader python class. Generate with,
readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')
Here is the reader,
read_file('trackingActionReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'turn', 'volumeUID', 'status', 'p', 'e', 'x', 'y', 'z', 'px', 'py', 'pz']
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for
# gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Only accept the primary muon death and its first generation daughters birth
if (e.trackID == 1 and e.status == 1) or (e.parentTrackID == 1 and e.status == 0):
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID,
e.parentTrackID, e.turn, e.volumeUID, e.status, e.p, e.e,
e.x, e.y, e.z, e.px, e.py, e.pz ])
return True
Make the reader class and object
trackingActionReaderC <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader
We have many files to process. Let’s do this reading in parallel.
We need to create a reader object per parallel job
trackingActionReaders <- map(1:length(arrXrdFiles), function(i) trackingActionReaderC(artInputTag("artg4")))
Load the data
jobs <- lapply(1:length(arrXrdFiles), function(i) getGalleryData(arrXrdFiles[i], trackingActionReaders[[i]]) %>%
mcparallel(name=i))
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Closing file, read 663158402 bytes in 161 transactions
Timings: Overall time = 200.338 s
Time per event: min=0.000 mean=0.019 max=7.364
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root
Closing file, read 657073380 bytes in 158 transactions
Timings: Overall time = 228.757 s
Time per event: min=0.000 mean=0.022 max=27.028
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root
Closing file, read 668231087 bytes in 160 transactions
Timings: Overall time = 208.320 s
Time per event: min=0.000 mean=0.020 max=16.990
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 664723747 bytes in 161 transactions
Timings: Overall time = 200.955 s
Time per event: min=0.000 mean=0.019 max=8.426
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root
Closing file, read 666316280 bytes in 160 transactions
Timings: Overall time = 234.668 s
Time per event: min=0.000 mean=0.022 max=17.103
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root
Closing file, read 664218464 bytes in 136 transactions
Timings: Overall time = 535.363 s
Time per event: min=0.000 mean=0.053 max=334.559
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root
Closing file, read 660527554 bytes in 161 transactions
Timings: Overall time = 266.787 s
Time per event: min=0.000 mean=0.026 max=21.746
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root
Closing file, read 656304716 bytes in 158 transactions
Timings: Overall time = 209.096 s
Time per event: min=0.000 mean=0.021 max=17.045
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root
Closing file, read 668176375 bytes in 160 transactions
Timings: Overall time = 209.131 s
Time per event: min=0.000 mean=0.020 max=11.414
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 664807569 bytes in 161 transactions
Timings: Overall time = 216.531 s
Time per event: min=0.000 mean=0.021 max=17.764
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root
Closing file, read 666316217 bytes in 160 transactions
Timings: Overall time = 227.132 s
Time per event: min=0.000 mean=0.022 max=13.152
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root
Closing file, read 665617262 bytes in 161 transactions
Timings: Overall time = 221.986 s
Time per event: min=0.000 mean=0.021 max=9.361
trackingActionReturns <- mccollect(jobs)
Let’s merge all of the output.
trackingActionDF <- map_df(1:length(trackingActionReaders),
function(i) galleryReader_df(trackingActionReaders[[i]]) %>% mutate(scenario=scenarios[i],
fileEntry=i))
trackingActionDF
write_rds(trackingActionDF, 'trackingActionDF.rds')
How many muons per file?
trackingActionDF %>% filter(trackID==1) %>% group_by(scenario, fileEntry) %>% tally()
We need to turn the volume IDs into volume names. The volume IDs change for each file. We can run a small art job to determine the volume ID to name tables.
runForVolIDString <- function(i) {
str_interp(
"PVS_CSVOUT=${csvout}_${i}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
list( csvout=scenarios[i], i=i, fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrXrdFiles[i])
)
}
runForVolIDStrings <- sapply(1:length(arrXrdFiles), runForVolIDString)
runForVolIDStrings
[1] "PVS_CSVOUT=everything_1_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"
[2] "PVS_CSVOUT=noDipole_2_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"
[3] "PVS_CSVOUT=noInflector_3_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"
[4] "PVS_CSVOUT=noKickerQuads_4_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
[5] "PVS_CSVOUT=noKicker_5_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"
[6] "PVS_CSVOUT=noQuads_6_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"
[7] "PVS_CSVOUT=everything_7_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"
[8] "PVS_CSVOUT=noDipole_8_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"
[9] "PVS_CSVOUT=noInflector_9_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"
[10] "PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "PVS_CSVOUT=noKicker_11_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"
[12] "PVS_CSVOUT=noQuads_12_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"
Let’s run the art jobs in parallel
jobs <- lapply(1:length(arrXrdFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)
running command 'PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root' had status 20
and load the CSV files.
jobs <- lapply(1:length(arrFiles), function(i)
str_interp("${scn}_${i}_volNames.csv", list(scn=scenarios[i], i=i)) %>%
read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
#
volNameTables <- mccollect(jobs)
volNameDF <- map_df(1:length(volNameTables), function(i)
volNameTables[[i]] %>% mutate(scenario=scenarios[i],
fileEntry=i))
volNameDF
write_rds(volNameDF, 'volNameDF.rds')
Now we merge with the tracking action data frame
trackingActionDF %>% inner_join(volNameDF) %>%
select(scenario, fileEntry, eventEntry, volName, everything()) -> trackingActionNDF
Joining, by = c("fileEntry", "volumeUID", "scenario")
trackingActionNDF
These muons decay in the storage region (not quite the same as “stored”),
trackingActionNDF %>% filter(trackID==1, volName=='ArcSection[00]') %>%
group_by(scenario, fileEntry) %>% tally()
Interesting! So, from this, only 0.575% decay in the storage region for the “all on” (or everything) sample. That seems rather low. But continuing.
write_rds(trackingActionNDF, 'trackingActionNDF.rds')
trackingActionNDF %>% filter(trackID == 1) %>% group_by(scenario, volName) %>% tally()
It may be nice to categorize the volume names by area (e.g. “Arc”, “bellows”, “Ring”, “Inflector”).
# Pull out the first word in camel case (see http://stackoverflow.com/questions/29916065/how-to-do-camelcase-split-in-python)
camelCaseSplit <- function(s) str_split(s, '(?<=[a-z])(?=[A-Z])|(?<=[A-Z])(?=[A-Z][a-z])')
# Test the above
expect_equal(camelCaseSplit("RingPoleTipUpper")[[1]], c("Ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("ringPoleTipUpper")[[1]], c("ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("world")[[1]], c("world"))
expect_equal(camelCaseSplit("ringPole")[[1]], c("ring", "Pole"))
# We really just want the first word
takeFirstCamel <- function(s) camelCaseSplit(s) %>% map_chr(function(s) s[1])
takeFirstCamel("RingPoleTipUpper") %>% expect_equal("Ring")
takeFirstCamel("ringPoleTipUpper") %>% expect_equal("ring")
takeFirstCamel("world") %>% expect_equal("world")
takeFirstCamel("ringPole") %>% expect_equal("ring")
takeFirstCamel(c("ringPoleTipUpper", "world")) %>% expect_equal(c("ring", "world"))
Don’t confuse RingTrackingPlane with the iron
makeVolCategory <- function(s) ifelse( str_detect(s, 'RingTracking'),
"RngTrkPlane",
takeFirstCamel(s))
makeVolCategory(c("RingBottom", "RingTrackingPlane[8]", "world")) %>% expect_equal(c("Ring", "RngTrkPlane", "world"))
trackingActionNDF %>% mutate(volCategory = volName %>% makeVolCategory()) -> trackingActionNVDF
trackingActionNVDF
Show count by category
trackingActionNVDF %>% filter(trackID == 1) %>% group_by(scenario, volCategory) %>% tally() -> volCatTally
volCatTally %>% mutate(perc = n/20000 * 100) -> volCatTally
volCatTally
Let’s plot!
trackingActionNVDF %>% filter(trackID ==1) %>% ggplot(aes(x=volCategory, group=scenario)) +
geom_bar()
ggplot(volCatTally, aes(x = volCategory, y=perc, fill=scenario)) +
geom_bar(stat="Identity", width=0.5, position="dodge")
Can’t really see non-Ring/world ones. Let’s split them out.
# Function to make this easy
plotit <- function(d) ggplot(d, aes(x = volCategory, y=perc, fill=scenario) ) +
geom_bar(stat="Identity", width=0.5, position="dodge") +
xlab("Volume category") + ylab("Percent")
volCatTally %>% filter(volCategory %in% c("Ring", "world")) %>% plotit -> p1
volCatTally %>% filter(! volCategory %in% c("Ring", "world")) %>% plotit -> p2
grid.arrange(p1, p2)
Let’s see where they die
# See https://rpubs.com/hadley/97970 for how to wrap a multipart ggplot2 component
plotCommon <- function () {
list(
facet_wrap(~scenario),
guides(col = guide_legend(override.aes = list(size=5, alpha=1))),
scale_color_discrete(name="Volume\nCategory"),
labs(x="z (mm)", y="x (mm)"),
theme_minimal()
)
}
trackingActionNVDF %>% filter(trackID == 1) %>%
ggplot( aes(x=z, y=x, color=volCategory)) +
geom_point(alpha=0.1) +
ggtitle('Where Muons Die') +
plotCommon()
Here it is without the ring losses, which dominate
trackingActionNVDF %>% filter(trackID == 1, volCategory != 'Ring') %>%
ggplot( aes(x=z, y=x, color=volCategory)) +
geom_point(alpha=0.5) +
ggtitle("Where Muons Die", subtitle = "Ring escapees excluded for clarity") +
plotCommon()
We have 3 GeV muons. Does it make sense that the vast majority of them stop in the iron? A table shows the Continuous Slow Down Approximation range of muons (CSDA). For a 3 GeV muon in iron, the CSDA is \(1.825 \times 10^{3}\) g/cm\(^2\). The density \(\rho\) is 7.874 g/cm\(^3\). So, where \(R\) is range,
\[R = \text{CSDA}/\rho\] Range is thus 2.32 m. This is begger than the width of the iron, but remember that the muons come in at an oblique angle. We’ll have to prove this.
We have the Geant stepping information, so we should be able to figure out how far the muons go in the iron.
Here is the Gallery reader,
read_file('GeantTrackReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class GeantStepReader(GalleryReaderBase):
def __init__(self, inputTag, eventsToFill = []):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'totalEnergyDeposit',
'stepLength', 'volumeUID',
'x', 'y', 'z', 'px', 'py', 'pz']
self.eventsToFill = eventsToFill
print self.eventsToFill
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))
def fill(self, ROOT, ev):
# Do we care about this event - if not then don't bother looking?
if self.eventsToFill and not ev.eventEntry() in self.eventsToFill:
return True
print 'Reading entry %s' % ev.eventEntry()
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2ringsim::GeantTrackRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2ringsim::GeantStep
for e in p: # Loop over elements and fill
# Get the steps
for f in e.steps():
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(),
f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(),
f.pos()[0], f.pos()[1], f.pos()[2],
f.p()[0], f.p()[1], f.p()[2] ])
# If we're on the last event, don't bother reading any more
if self.eventsToFill and ev.eventEntry() == self.eventsToFill[-1]:
return False
return True
geantStepReaderC <- createReaderClass_from_file('GeantTrackReader.py')$GeantStepReader
Let’s find an event where the muon did in the iron.
trackingActionNVDF %>% filter(trackID == 1, volName == 'RingYokeBottom')
Event 4 looks good.
geantStepOne <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(4)))
(4,)
getGalleryData(arrXrdFiles[1], geantStepOne)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 4
gsOneDF <- galleryReader_df(geantStepOne)
Merge in volume names
gsOneDF %>% mutate(fileEntry = 1) %>%
inner_join(volNameDF) %>%
select(volName, everything()) -> gsOneDF
Joining, by = c("fileEntry", "volumeUID")
gsOneDF
Let’s make a 3D plot of this path
ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
radius=7112,
sides=20, closed = F)
Let’s do volume categories again.
gsOneDF %>% mutate(volCategory = volName %>% makeVolCategory %>% as.factor) -> gsOneDF
Let’s plot the path. Note that I have to make the legend separately since legend3d looks awful.
# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneDF,
legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsOneDF,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
How much energy loss and distance traveled in each volume category
gsOneDF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n())
Ok - I believe this (note that the y scale is very small compared to the x,z scales - so the muon acutally travels quite far in x & z). The muon travels about 2.4m in the iron and ranges out. Physics (and Geant) work!
Let’s try another one
geantStepOneA <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(10)))
(10,)
getGalleryData(arrXrdFiles[1], geantStepOneA)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 10
Timings: Overall time = 178.569 s
Time per event: min=0.000 mean=0.100 max=1.101
gsOneADF <- galleryReader_df(geantStepOneA)
gsOneADF %>% mutate(fileEntry = 1) %>%
inner_join(volNameDF) %>%
mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>%
select(volName, volCategory, everything()) -> gsOneADF
Joining, by = c("fileEntry", "volumeUID")
gsOneADF
gsOneADF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsOneADFSum
gsOneADFSum
Check sums
gsOneADFSum %>% summarize(sum(totalELoss_MeV))
gsOneADF %>% summarize(sum(totalEnergyDeposit))
# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneADF,
legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsOneADF,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
This one scatters in the inflector.
Can we show how far muons go in various materials?
Let’s collect data for 100 events.
We have a reader that can do n events.
read_file('GeantTrackReaderN.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class GeantStepReaderN(GalleryReaderBase):
def __init__(self, inputTag, startEvent=100, nEvents=10):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'globalTime', 'totalEnergyDeposit',
'stepLength', 'volumeUID',
'x', 'y', 'z', 'px', 'py', 'pz']
self.startEvent = startEvent
self.endEvent = startEvent + nEvents
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))
def fill(self, ROOT, ev):
# Do we care about this event - if not then don't bother looking
if ev.eventEntry() < self.startEvent:
return True
if ev.eventEntry() >= self.endEvent:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2ringsim::GeantTrackRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2ringsim::GeantStep
for e in p: # Loop over elements and fill
# Get the steps
for f in e.steps():
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), f.globalTime(),
f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(),
f.pos()[0], f.pos()[1], f.pos()[2],
f.p()[0], f.p()[1], f.p()[2] ])
return True
geantStepReaderNC <- createReaderClass_from_file('GeantTrackReaderN.py')$GeantStepReaderN
geantStepN <- geantStepReaderNC(artInputTag('artg4:KeepStepsAction'), 1000, 1000)
getGalleryData(arrXrdFiles[1], geantStepN)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Timings: Overall time = 66.083 s
Time per event: min=0.000 mean=0.030 max=2.367
gsNDF <- galleryReader_df(geantStepN)
gsNDF %>% mutate(fileEntry = 1,
p = sqrt(px*px + py*py + pz*pz)) %>%
inner_join(volNameDF) %>%
mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>%
select(volName, volCategory, p, everything()) -> gsNDF
Joining, by = c("fileEntry", "volumeUID")
gsNDF
Let’s do a distribution of distance in the different volumes
gsNDF %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDFSum
gsNDFSum
gsNDFSum %>% filter(volCategory == 'Ring') %>%
ggplot(aes(x=totalDist_meters)) + geom_histogram(bins=50) +
labs(x = 'Total distance in Ring material (m)', y='count', title='Muon in Ring Material')
Let’s look at this across the volume categories for energy loss
gsNDFSum %>% ggplot(aes(x=totalELoss_MeV)) + geom_histogram(bins=50) +
facet_wrap(~volCategory, scales = "free") +
labs(x = 'Total energy loss (MeV)', y='count', title='Energy loss in material',
subtitle = 'Note the different scales')
Let’s get an energy loss curve
gsNDF %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) +
labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))")
Well, not sure what all this means. Let’s look for a particular event…
gsNDF %>% filter(eventEntry == 1300) %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) +
labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))", title='Energy loss vs. Momentum for muons in event 1300' )
How many steps do we take?
gsNDFSum %>% ggplot(aes(x=nSteps)) + geom_histogram(bins=50) +
facet_wrap(~volCategory, scales = "free") +
labs(x = 'Number of steps', y='count', title='Number of steps in materials',
subtitle = 'Note the different scales')
Let’s look at two events
gsNDFSum %>% filter(volCategory == 'Ring', totalDist_meters > 3)
Let’s look at event 1758
gsNDF %>% filter(eventEntry == 1758) -> gsNDF1758
gsNDF1758
gsNDF1758 %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDF1758Sum
gsNDF1758Sum
Let’s plot it!
gsNDF1758 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1758
# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1758,
legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsNDF1758,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
Looks like this one started out in a strange place.
gsNDFSum %>% filter(volCategory == 'Arc', totalDist_meters > 870)
Let’s try event 1509
gsNDF %>% filter(eventEntry == 1509) -> gsNDF1509
gsNDF1509
gsNDF1509 %>% group_by(eventEntry, volCategory) %>% summarize(
totalDist_meters=round(sum(stepLength)/1000, 2),
totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
nSteps = n()) -> gsNDF1509Sum
gsNDF1509Sum
Let’s plot it!
gsNDF1509 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1509
# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1509,
legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)
clear3d()
with(gsNDF1509,{
plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
Looks like a nice stored muon!! Can we see CBO?
gsNDF1509 %>% mutate(r = sqrt(x*x+z*z)) -> gsNDF1509
p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + geom_line() +
labs(x='Global time (ns)', y='Radius (mm)', title='Radial CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) +
geom_line() +
labs(x='Global time (ns)', y='Radius (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)
We can even see the kick!
Here we’ll plot the points too.
p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + geom_line() +
geom_point() +
labs(x='Global time (ns)', y='y (mm)', title='Vertical CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) +
geom_line() + geom_point() +
labs(x='Global time (ns)', y='y (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)